Application  of  MCT  failure  criterion 

USING  EFM 


1  Summary 

This  document  describes  the  work  done  in  NUS  to  implement  the  Multi-Continuum 
Theory  (MCT)  of  composite  failure  for  damage  progression  analysis  using  the  Element-Failure 
Method  (EFM)  in  ABAQUS.  The  EFM-MCT  has  been  implemented  in  ABAQUS  using  a  user 
defined  element  subroutine  EFM.  Comparisons  have  been  made  between  the  analysis  results 
using  EFM-MCT  code  and  HELIUS:MCT™  code  which  uses  material  property  degradation 
method  (MPDM).  It  is  found  that  the  conventional  EFM-MCT  code  is  less  efficient  than  the 
HELIUS:MCT™  code.  This  is  due  to  the  iterative  nature  of  the  EFM  code  in  finding  the 
solution.  A  non-iterative  element-failure  method  (N-EFM)  is  then  developed  to  improve  the 
computational  efficiency.  This  N-EFM-MCT  has  been  implemented  in  an  in-house  code  and 
comparison  of  computational  efficiency  has  been  made  between  N-EFM  and  MPDM. 

2  Multi-continuum  theory 

Multi-continuum  theory  (MCT)  [1-15]  is  a  micromechanics-based  failure  theory 
applicable  to  fiber  reinforced  plastic  which  is  based  on  phase  averaging  of  stress  and  strain. 
Based  on  the  averaging  of  stress  and  strain  of  constituent  matrix  and  fiber  in  a  composite,  a 
relationship  between  stress  and  strain  of  the  composite  and  the  stress  and  strain  of  the 
constituents  can  be  established,  i.e.  constituent  (matrix  and  fiber)  stress  and  strain  can  be 
calculated  from  the  composite  stress  and  strain.  The  failure  can  then  be  evaluated  based  on  the 
constituent  stress. 

The  phase  averaging  of  stress  and  strain  in  MCT  for  a  limited  representative  volume  is 


expressed  as  [1]: 

(1) 

(2) 
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where  and  ff^’^are  6x1  vectors  that  represent  respectively  composite  ,  fiber  and  matrix 

average  stress;  and  are  6  x  1  vectore  that  represents  composite  ,  fiber  and  matrix 

average  stress  matrices,  respectively;  while  and  are  respectively  fiber  and  matrix 
volume  fraction. 

MCT  uses  the  following  linear  elastic  constitutive  relationship  between  stress  and  strain: 


(3) 


(4) 


(5) 


where  C^,  Cr,  andCjjj^  represent  6x6  stiffness  matrices  of  composite,  fiber  and  matrix, 

respectively;  6  is  temperature  change;  cfp,  Cf^,  are  6  x  1  vectors  of  composite,  fiber  and 

matrix  thermal  expansion  coefficients,  respectively.  By  using  equations  (1)  to  (5),  relationships 
between  composite  stress  and  strains  with  constituent  (fiber  and  matrix)  stress  and  strains  can  be 
determined;  hence  constituent  stress  and  strain  can  be  calculated  from  composite  stress  and 
strain. 

Once  the  constituent  stress  and  strain  are  determined,  MCT  evaluates  composite  failure 
based  on  its  constituent  stress.  The  failure  criterion  used  in  MCT  for  matrix  can  be  expressed  as 


[1]: 


(6) 


while  the  failure  criterion  for  fiber  is 


(V) 


where 


(8) 


(9) 
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/r= 


(10) 


(11) 


(12) 


(8) 


and  A^  are  failure  coefficients  that  can  be  derived 
from  tensile,  compression,  and  shear  test  data  of  the  composite  [1], 

3  Element  failure  method 

This  section  describes  the  conventional  iterative  element-failure  method  (EFM)  and  the 
modified  non-iterative  element-failure  method  (N-EFM). 

3.1  Iterative  Element-Failure  Method 

Element-failure  method  (EFM)  is  a  damage  progression  algorithm  in  which  the  effect  of 
damage  in  mechanical  behavior  is  affected  through  nodal  forces  in  finite  element  (FE)  analysis. 
The  concept  was  first  conceived  for  dynamic  fracture  in  metals  [16],  but  the  modified  EFM  was 
found  useful  for  the  analysis  of  impact  damage  in  fiber-reinforced  composites  [17],  damage 
progression  in  quasi-statically  loaded  three  point  bend  composite  laminates  [18]  and  ultimate 
strengths  of  open-hole  tension  composite  laminates  [19].  The  manner  by  which  these  effects  due 
to  damage  translate  to  the  effective  nodal  forces  will  in  general  depend  upon  the  damage 
evolution  law  appropriate  to  the  local  mode  of  damage  experienced  by  the  composite  material,  as 
well  as  the  FE  formulation.  Traditionally,  when  damage  is  assumed  to  have  occurred  within  an 
element  of  a  material,  it  is  reasonable  to  assume  that  the  stiffness  matrix  of  the  element  is  altered 
to  reflect  the  damaged  state.  As  the  following  derivation  shows,  there  are  explicit  relations 
between  the  nodal  forces  and  the  elastic  stiffness  of  an  FE. 

The  force-stiffness  relation  for  an  FE  is  given  by 


Klf  =  f 


(9) 
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where  u  is  the  vector  of  nodal  displacements,  /  the  vector  of  nodal  forces  and  the  element 
stiffness  matrix  of  undamaged  material.  In  FE  analysis,  the  stiffness  matrix  K  is  calculated  by 
integrating  over  the  element  domain: 

P  (10) 

if  =  I  B^CB  dn 

where  B  is  the  derivative  of  element  shape  function  and  C  is  the  material  stiffness  matrix. 

In  conventional  progressive  damage  algorithm  which  uses  material  property  degradation 
method  (MPDM),  element  or  material  failure  is  modeled  by  modifying  the  material  stiffness  C 

and  hence  the  element  stiffness  matrix  K.  EFM,  however,  models  the  failure  by  directly 

modifying  the  force  vector  /.  Consider  an  FE  of  an  undamaged  composite  material  (Figure  1 

(a)),  experiencing  a  set  of  nodal  forces,  which  have  been  obtained  from  the  FE  solution  of  the 
problem.  On  the  other  hand,  in  an  FE  containing  damaged  material,  the  load-bearing  capacity  of 
the  FE  will  be  compromised,  very  likely  in  a  directionally  and  spatially  dependent  manner.  If 
much  of  the  damage  consists  of  transverse  matrix  microcracks,  it  is  reasonable  to  assume  that  the 
FE  will  have  reduced  load-bearing  capacity  in  the  direction  transverse  to  the  fibers  (Figure  1  (b)). 
In  conventional  MPDM,  this  reduction  is  achieved  by  reducing  or  zeroing  certain  pertinent 
material  stiffness  properties  of  the  damaged  FE.  In  the  EFM,  however,  the  reduction  is  effected 
by  applying  a  set  of  external  nodal  forces  such  that  the  net  internal  nodal  forces  of  elements 
adjacent  to  the  damaged  element  are  reduced  or  zeroed  (the  latter  if  complete  failure  or  fracture 
is  implied  (Figure  1  (c)).  The  decision  whether  to  fail  an  element  is  guided  by  a  suitable  failure 
theory.  The  required  set  of  applied  nodal  forces  to  achieve  the  reduction  within  each  step  is 
determined  by  successive  iterations  until  the  net  internal  nodal  forces  (residuals)  of  the  adjacent 
elements  converge  to  the  desired  values.  Note  that  it  is  not  the  internal  nodal  forces  of  the 
damaged  element  that  is  zeroed  (for  the  case  of  complete  failure  (Figure  1  (c)),  but  the  net 
internal  nodal  forces  of  adjacent  elements.  Thus  the  ‘stresses’  within  the  failed  element  no  longer 
have  physical  meaning,  although  compatibility  may  be  preserved.  This  process  leaves  the 
original  (undamaged)  material  stiffness  properties  unchanged,  and  is  thus  computationally 
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efficient  as  every  step  and  iteration  are  simply  analysis  with  the  updated  set  of  applied  nodal 
forces.  For  this  reason,  it  may  also  be  called  the  nodal  force  modification  method.  Hence,  no 
reformulation  of  the  FE  stiffness  matrix  is  necessary. 


Figure  1  (a)  FE  of  undamaged  composite  with  internal  nodal  forces  ;  (b)  FE  of  composite  with 
transverse  matrix  cracks.  Components  of  internal  nodal  forces  transverse  to  fiber  direction  are 
modified;  (c)  Completely  failed  or  fractured  element.  All  net  internal  nodal  forces  of  adjacent 

elements  are  zeroed 

3.2  Non-Iterative  Element-Failure  Method  (N-EFM) 

Non-iterative  element  failure  method  (N-EFM)  is  focused  on  implementation  of  damaged 
elements  with  direct  solution  and  implicit  external  nodal  forces  applied  to  the  damaged  elements. 
Compared  with  the  traditional  FEM  with  MPDM,  N-EFM  doesn’t  require  frequent  assembly  and 
decomposition  of  global  stiffness  matrix  in  each  load  increment.  The  basic  equation  in  N-EFM  is 
derived  from  the  general  FEM  as: 


(16a) 


Kuj  =  F, 


where 


(17) 


and  Kq  is  the  initial  stiffness  matrix  without  damage,  which  keeps  constant  through  the  analysis; 
AK  is  the  stiffness  matrix  for  damaged  elements,  which  is  updated  with  failure  progression;  Uj 
and  Pj  are  the  displacements  and  known  loads  in  the  /th  increment. 

In  Eq.(17),  only  one-time  LU  decomposition  (e.g.  Cholesky  factorization,  =  LlJ)  at 


the  beginning  of  the  analysis  is  required  and  then  stored  for  use  in  subsequent  computations.  Fj 
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can  be  easily  obtained  by  forward  and  backward  substitutions  with  the  LU  decomposition. 
Additionally,  is  also  not  explicitly  calculated  and  only  columns  related  to  damaged  elements 

are  required,  which  are  obtained  by  a  same  method  with  that  of  obtaining  Fj .  The  multiplication 
fir“'AK'  in  Eq.(17)  accordingly  takes  the  form  as: 


0  • 

Kj 

■  0“ 

0  • 

••  •• 

•  0 

(18) 

0  • 

Kj 

■  0 

(nxn) 

where  only  columns  related  to  damaged  elements  are  non-zero. 

As  is  constant  through  the  analysis,  only  columns  in  related  to  new  damaged 
elements  are  calculated  in  each  load  increment  and  then  stored  for  use  in  subsequent  load 
increments.  After  K  and  Fj  are  obtained,  Eq.(16a)  can  be  rearranged  with  the  following 
formulation: 


1 _ 

o’ 

1 

1 _ 

< 

1 _ 

_ 1 

1 

1 - 

s 

1 _ 

1 

_ 1 

(16b) 


where  and  are  condensed  matrices  related  to  damaged  and  undamaged  elements, 

respectively,  and  the  subscript  ’  and  “w”  denote  damaged  elements  and  undamaged  elements, 
respectively. 

Then  the  large-size  problem  described  by  Eq.(16a)  or  Eq.(16b)  can  be  reduced  to  a  small- 
size  problem  described  by  the  following  equation: 

(19) 

where  the  solution  is  only  focused  on  damaged  elements  and  the  problem  size  is  dependent  on 
the  number  of  damaged  elements. 
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After  the  displacement  related  to  damaged  elements,  ,  is  obtained,  the  displacement 
related  to  undamaged  elements,  ,  can  be  obtained  straightforward  with  substituting  into  the 
following  equation  according  to  Eq.(16b): 

(20) 

In  N-EFM,  the  global  stiffness  matrix  for  damaged  elements,  AK ,  as  indicated  in 
Eq.(15),  is  assembled  with  the  element  stiffness  matrix  for  damaged  elements,  AK\  which  is 
determined  by  the  following  equations  with  MCT: 

\k%  if  fibre  fails 

AK^=\  (21) 

[K^ ,  if  matrix  fails 


=  yB^D  B 

m  m 


(22) 


where  D  and  are  full  and  reduced  global  material  stiffness  matrices,  respectively;  D  and 
are  obtained  with  transformation  of  full  and  reduced  local  material  stiffness  matrices,  D  and 
,  respectively,  and 


jo,  if/,7=l 

[Dy, 


(23) 


In  order  to  speed  up  the  computation  of  N-EFM,  the  initial  stiffness  matrix  without 
damage  in  Eq.(17),  K^,  can  be  replaced  by  a  reassembled  stiffness  matrix  after  every 

increments,  Kj^ ,  and  AK  is  also  replaced  by  that  after  every  increments,  AKj^ ,  which  is 
assembled  only  for  new  damaged  elements  in  every  increments.  Then  the  N-EFM  equation 
can  be  modified  as: 

(24) 


and  the  reassembled  stiffness  matrix,  Kj^ ,  is  determined  by  the  following  equation: 
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(25) 


K, 


nit  ~^d 

'Z.K-+Y.K- 


where  m,  and  are  total  number  of  elements  and  number  of  damaged  elements,  respectively, 
and 

Ik'%  if  fibre  fails 

K=\  (26) 

[fSr;;;,  if  matrix  fails 


and  the  element  stiffness  matrix  K'‘'  and  K'^  are  determined  by  the  local  material  stiffness 
matrices  Z)'.  and  with  appropriate  transformations,  which  are  defined  as: 


(27) 


jZ)..,  ifi,7  =  l 

[dj^Dy, 


(28) 


where  w  0  is  a  material  degradation  factor. 


4  Implementation  of  EFM-MCT  in  ABAQUS 

EFM-MCT  has  been  implemented  in  ABAQUS  using  user  defined  element  subroutine 
UEL.  The  exactly  same  MCT  failure  criterion  used  in  the  HELIUS:MCT™  code  [1]  is  used  in 
the  EFM-MCT  code.  In  this  code,  failure  of  each  element  is  examined  in  its  centroid.  EFM  nodal 
force  modification  is  performed  whenever  any  element  is  considered  failed.  If  matrix  in  an 
element  is  considered  failed,  the  nodal  force  is  modified  so  that  nodal  force  due  to  E2,  E'3,  En, 
Eu,  and  E23  are  negated.  If  fiber  in  an  element  is  considered  failed,  the  nodal  force  is  modify 
such  that  it  negates  the  stiffness  of  the  element. 

Damage  analysis  using  MCT  requires  the  calculation  of  matrix,  fiber,  and  composite 
coefficients  which  include  stiffness,  Poisson’s  ratio,  and  failure  coefficients  before  the  analysis 
begins.  HELIUS:MCT™  uses  Helius  material  manager  to  calculate  these  coefficients  based  on 
the  measured  composite  stiffness,  Poisson’s  ratio  and  strength,  and  the  measured  constituent 
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(matrix  and  fiber)  stiffness  and  Poisson’s  ratio.  The  EFM-MCT  code  that  is  developed  in  this 
study  also  uses  Helius  material  manager  to  calculate  the  coefficients. 

Damage  analysis  using  EFM-MCT  requires  two  simulations:  the  first  simulation  is  the 
actual  progressive  damage  analysis,  while  the  second  simulation  is  dummy  analysis  for 
visualization  purpose.  The  second  simulation  is  needed  because  ABAQUS  does  not  support 
visualization  for  elements  from  UEL  subroutine.  Thus,  in  this  code,  element  failure  information 
from  the  E'  analysis  is  stored  in  a  text  file.  This  text  file  is  then  read  in  the  second  analysis  and 
the  element  failure  information  is  stored  in  ABAQUS  native  elements  so  that  it  can  be  visualized 
in  ABAQUS. 

5  Comparison  between  HELlUSiMCT™  and  EFM-MCT 

Open  hole  tension  (OHT)  models  are  analyzed  using  HELIUS:MCT™  and  EFM-MCT 
codes  in  ABAQUS  in  order  to  compare  the  two  methods.  The  plate  is  76.2  mm  x  76.2  mm  in  size 
and  it  has  a  19.05  diameter  central  hole.  Carbon/epoxy  IM7/5250-4  composite  laminate  with 
stacking  sequence  [45/0/-45/90]s  is  used  for  the  plate  and  their  properties  are  listed  in  Table  1. 
Four  different  finite  element  meshes  were  used  to  study  the  sensitivity  of  these  codes  to  mesh 
size.  Figure  2  shows  four  different  finite  element  meshes  with  168,  400,  600  and  864  elements, 
respectively. 

Thermal  analysis  is  applied  in  the  simulation  using  EFM-MCT  code  in  order  to  take  into 
account  the  ply  thermal  residual  stress  that  occurs  during  the  manufacturing/curing  process  of  the 
composite.  Simulations  using  HELIUS  :MCT™  on  the  other  hand  do  not  include  thermal 
analysis  because  HELIUS:MCT™  does  not  facilitate  this.  Attempts  have  been  made  to  use 
ABAQUS  native  thermal  expansion  model  combined  in  addition  to  Helius-MCT  material  model 
in  order  to  model  the  curing  process,  but  a  converged  solution  cannot  be  obtained.  Thus,  another 
simulation  using  EFM-MCT  without  the  thermal  analysis  is  performed  so  that  the  Helius-MCT 
can  be  compared  to  EFM-MCT  with  the  exactly  same  parameters. 

The  simulation  results  in  the  average  stress  strain  curves  are  shown  in  Figure  3.  They 
suggest  that  both  EFM-MCT  and  HELIUS:MCT™  codes  are  not  sensitive  to  mesh  size  when 
the  mesh  is  reasonably  fine,  i.e.  more  than  400  elements  in  this  case.  Moreover,  comparison 
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between  predicted  ultimate  strength  with  experimental  result  [20],  as  shown  in  Figure  4,  shows 
that  both  models  are  able  to  predict  the  experiment  tensile  strength  quite  well  with  error  less  than 
15%.  Simulations  using  EFM-MCT  that  do  not  include  thermal  analysis  predict  a  slightly  higher 
strength  than  the  results  predicted  by  simulations  that  include  thermal  analysis;  their  difference  is 
around  1.3%  -  3.3%. 


Table  1  Material  properties  of  IM7/5250-4  composite  system  [20] 


Thickness  (mm) 

0.125 

Modulus  in  fiber  direction  E\  (GPa) 

172.4 

Transverse  moduli  £2=  E3  (GPa) 

10.3 

Shear  moduli  G\2=Gi3  (GPa) 

5.52 

Shear  modulus  G23  (GPa) 

3.45 

Poisson’s  ratios  Vi2=  vb 

0.32 

Poisson’s  ratio  V23 

0.4 

Longitudinal  tensile  strength  X  (MPa) 

2826.5 

Longitudinal  compressive  strength  X’  (MPa) 

1620.0 

Transverse  tensile  strength  Y  (MPa) 

65.5 

Transverse  compressive  strength  Y’  (MPa) 

248.0 

Shear  strength  5'i2=*S'i3  (MPa) 

122.0 

Shear  strength  S23  (MPa) 

122.0 

Coefficient  of  thermal  expansion  in  fiber  direction  ai(/°F) 

-0.2x10'^ 

Coefficients  of  thermal  expansion  in  transverse  directions  02=  03 
(/°F) 

13.7x10'*^ 

Temperature  difference  (°F) 

-375 

Figure  2  Finite  element  mesh  for  open-hole  tension  model 
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HeliusMCT 


Strain 


Figure  3  Stress-strain  eurves 


Figures  4,  5  and  6  show  the  failure  pattern  when  major  load  drops  oecur  in  the 
simulations  using  EFM-MCT,  EFM-MCT  -  no  thermal,  and  HELIUS:MCTT'^,  respeetively. 
Figures  4  and  5  suggest  that  although  strengths  predieted  by  simulations  with  and  without 
thermal  analysis  are  almost  the  same,  their  failure  patterns  are  largely  different.  Simulations 
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using  thermal  analysis  predict  a  much  larger  matrix  failed  area  than  simulations  without  thermal 
analysis.  This  shows  that  thermal  residual  stress  accounted  for  in  the  simulations  using  thermal 
has  a  major  influence  on  matrix  failure. 


Figure  5  OHT  failure  pattern  using  EFM-MCT  with  thermal/curing  analysis 
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Figure  6  OHT  failure  pattern  using  EFM-MCT  without  thermal/curing  analysis 
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Figure  7  OHT  failure  pattern  using  HELIUS:MCTT'^  without  thermal/curing  analysis 


Analysis  using  Mesh  4  with  different  number  of  time  step  was  performed  to  compare  the 
computation  time  needed  by  HELIUS:MCT™  and  EFM  MCT.  Table  2  shows  the  comparison  of 
CPU  time  needed  by  the  two  codes.  Note  that  the  simulation  time  for  EFM-MCT  consists  of  two 
parts:  the  first  one  is  from  the  actual  analysis  while  the  second  is  from  the  dummy  simulation 
used  for  visualization  purpose.  Table  2  suggests  that  HELIUS:MCT™  is  more  efficient  in  terms 
of  CPU  time  compared  to  EFM-MCT. 
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Table  2  CPU  time  for  864  element  model 


HELIUS:MCTTM 

EFM-MCT 

101  time  step 

401  time  step 

101  time  step 

401  time  step 

U‘run  (actual 
simulation) 

691.72 

1293.4 

3987.4 

6910.8 

2“'^  run  (post 
processing  for 
visualization) 

N/A 

N/A 

73.2 

180.67 

total 

691.72 

1293.4 

4060.6 

7091.47 

6  Comparison  between  N-EFM  and  MPDM 

An  example  of  crack  propagation  is  used  to  demonstrate  the  computational  efficiency  of 
N-EFM,  which  is  also  compared  with  that  of  MPDM.  The  model  and  mesh  are  shown  in  Figure 
8,  and  total  3690  3D  linear  brick  elements  with  20500  degree  of  freedoms  (DOFs)  are  used  in  the 
analysis. 


M  t  t  t  t  t 

E, 


-4- 


w 


I  I  felled  eleineubi 
I  I  mifeiled  elements 


Figure  8  Example  of  crack  propagation:  (a)  model;  (b)  mesh 


Figure  9  shows  the  comparison  of  compuational  efficiency  of  N-EFM  and  MPDM,  with 
normalized  compuation  time  by  that  of  MPDM,  where  N-EFM_1  and  N-EFM_2  denote  results 
from  reassembling  and  decomposing  global  stiffness  matrix  every  5  and  10  load  increments. 
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respectively,  and  3  elements  fail  in  each  load  increment.  It  can  be  seen  that  the  efficiency  of  N- 
EFM  is  improved  up  to  about  70%  compared  to  that  of  MPDM. 


- MPDM  -^N-EFM  1  ^leN-EFM  2 


Figure  9  Comparison  of  compuational  efficiency  of  N-EFM  and  MPDM 

7  Conclusion 

MCT  failure  criterion  has  been  implemented  in  ABAQUS  using  the  iterative  element- 
failure  method  (EFM)  algorithm.  OHT  analysis  is  performed  using  the  newly  developed  EFM- 
MCT  and  HELIUSrMCT^'^  which  uses  material  property  degradation  (MPDM)  algorithm. 
Results  show  that  EFM-MCT  and  HELIUS:MCT™  give  similar  results  and  they  are  in  good 
agreement  with  experimental  result.  The  EFM-MCT  is  however  less  efficient  than  the  MPDM 
HELIEFS:MCT^'^.  The  non-iterative  element- failure  method  (N-EFM),  which  is  a  modified 
version  of  conventional  iterative  EFM,  shows  better  computational  efficiency  compared  to 
standard  MPDM  when  the  number  of  failed  elements  are  small  in  large  finite  element  models. 
When  the  stiffness  matrix  is  reassembled  after  several  load  increments,  N-EFM  can  even  speed 
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up  more  computational  efficiency  (improved  up  to  70%).  Demonstrated  concept  so  far  is 
available  for  crack  problems.  We  have  coded  MCT  in  our  N-EFM  progressive  failure  code.  More 
work  on  demonstration  of  improvements  in  efficiency  with  OHT  problems  is  on-going. 
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